#####################
#MAP RUSSIAN TSARDOM#
#####################
#This to obtain:
#Figure_A20.jpg

rm(list = ls())
library("ggplot2")
library("dplyr")
library("ggrepel")
library("RColorBrewer")
library("lemon")
library("maptools")
library("ggsn")
library("ggplot2")
library("dplyr")
library("lemon")
library("raster")
library("readxl")
library("sf")

setwd("/Users/bgpopescu/Dropbox/")
#setwd("C:/Users/bogdanp/Dropbox/")

#Reading polygons
sf::sf_use_s2(FALSE)
all_coun <- st_read(dsn="./Legacies_Project/Analysis/data/data.gdb",
                    layer="all_countries_GCS")
all_coun<-st_simplify(all_coun,  dTolerance = 0.005)


colonies <- st_read(dsn="./Legacies_Project/Analysis/data/data.gdb",
                    layer="colonies")
colonies$y_lat<-sf::st_coordinates(st_centroid(colonies))[,2]
colonies$x_lon<-sf::st_coordinates(st_centroid(colonies))[,1]

#Reading rasters
euro_raster <- raster("./Legacies_Project/Analysis/data/euro_raster.tif")

#Preparing raster
max_lat<-max(colonies$y_lat) + 5
min_lat<-min(colonies$y_lat) - 5

max_lon<-max(colonies$x_lon) + 5
min_lon<-min(colonies$x_lon) - 5

e <- extent(min_lon-8, max_lon+8, min_lat-8, max_lat+8)
rc <- crop(euro_raster, e)	

#Lower resolution raster
rc.aggregate <- aggregate(rc, fact=3)

#Get rid of really high mountains
rc <- reclassify(rc.aggregate, c(-10000,0,NA))

#Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(rc); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Elev")


map<-ggplot()+
  geom_raster(data=hdf,aes(X,Y,alpha=Elev)) +
  scale_alpha(name = "Altitude", guide = "none")  + 
  geom_sf(data = all_coun, fill=NA, color = "black", size = 0.4)+
  geom_sf(data = colonies, fill="red", color = "red", size = 0.4)+
  coord_sf(ylim=c(min_lat,max_lat), xlim=c(min_lon, max_lon))+
  theme(legend.position="left")+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12),
        aspect.ratio = 1.2)+
  ggsn::scalebar(x.min = min_lon, x.max = max_lon,
                 #Above you mention how long the scalebar should be
                 y.min = min_lat+10, y.max = max_lat-2,
                 #Above you mention how thick the scalebar should be
                 dist = 100, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.05,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.05)

map2<-map+geom_label_repel(data=colonies, aes(label=name, x = x_lon, y = y_lat,
                                 fontface = 'bold'), alpha = 0.8)

ggsave(map2, file="./Legacies_Project/Paper/figures/Figure_A20.jpg", 
       height=20.2, width=20.2, units = "cm", dpi=300)


